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A binary mixture of interacting Bose-Einstein Condensates (BEC) in the presence of 
fragmentation-driving external lattice potentials forms two interdependent mean-field lattices made 
of each component. These effective mean-field lattices, like ordinary optical lattices, can induce 
fragmentation and phase coherence loss of BECs between lattice sites. We consider the nonequi- 
librium dynamics of two hyperfine states of ID Bose-Einstein condensates confined by a harmonic 
trap, subjected to state-dependent optical lattices. Our numerical mean-field calculations using the 
truncated Wigner approximation (TWA), include quantum and thermal fluctuations. We show loss 
of phase coherence in an unequal mixture of two-component BECs in a gradual ramp-up of one or 
two applied optical lattices, and we study coherence loss from mean-field lattices as a function of 
lattice height. 
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I. INTRODUCTION 

The properties of ultracold atoms in optical lattices 
have been studied intensively as models for various con- 
densed matter phenomena [l|, [|J and also as interest- 
ing quantum mechanical phenomena including quantum 
information and quantum computation [3- 6] . In many 
cases, the coherence or lack of coherence between atoms 
in adjacent wells plays a crucial role, especially in the 
superfluid-Mott insulator phase transition |7Hl0j. For 
example, the work of Orzel et al. [Tlj . with a ID array 
of "pancake" condensates, displayed high- visibility inter- 
ference patterns under conditions of phase coherence be- 
tween adjacent wells, but dramatic reduction of the inter- 
ference contrast, or visibility, when the wells were deep- 
ened and number-squeezing was believed to occur. Later 
work with 3D optical lattices by Greiner et al. ex- 
hibited the superfluid to Mott-Insulator phase transition 
through the interference pattern when the 3D conden- 
sates were released. 

For such single-component BECs in optical lattices, 
theoretical modeling and analysis of phase decoher- 
ence and number-squeezing phenomena in ID have been 
extensively developed both from an extended Gross- 
Pitaevskii equation (GPE) approach [ljj and from the 
Truncated Wigner Approximation (TWA) approach [l| . 
The application of the GPE to such systems is limited to 
shallow lattices and low temperatures unless used in a 
full 3D treatment that includes transverse excited modes 
[12j . The TWA, which evolved from quantum optics ap- 
plications 0, HU , has been subject to careful tests [l6[ , 
and has emerged as the method of choice for dealing with 
time evolution of Bose-Einstein condensates in optical 
lattices. The TWA has been used to model dephasing of 
(single component) BECs in ID optical lattices [IH, Il7l . 
Higher order corrections have been discussed in Ref. [l8[ 
but will not be considered here. 

More recently, experiments have addressed the ques- 
tion of phase coherence of systems in which the Bose con- 



densate consists of distinct components, such as different 
atomic species, different substates or different hyperfine 
levels. For example, it has been observed [19( that the 
presence of 41 K atoms reduces the visibility of the inter- 
ference pattern of marginally-overlapped 87 Rb atoms in 
a 3D optical lattice. Similarly, in a condensate of misci- 
ble 87 Rb atoms in a state-dependent 3D optical lattice, 
the presence of atoms in a second hyperfine level can re- 
duce the superfluid coherence of atoms in a first hyperfine 
level [13] ■ Lasers in these studies were tuned such that 
both components undergo peak-matched lattice poten- 
tials. These experiments are most closely related to the 
theoretical discussion in this study, but there are other 
recent experimental studies of Bose condensates focu sing 
on diverse aspects of multi-component systems [2ll425| 
and it is easy to imagine further experimental studies 
probing the rich physics of two-component BECs. 

Recently, also, there have been a number of theoretical 
studies addressing aspects of Bose condensates with two 
or more components. Perhaps most notable have been 
discussions of different phase regimes and phase transi- 
tions [26l - |28j and of the extended Bloch band structure 
[29l | . There have also been studies of dynamical effects 
including those associated with ramping up the optical 
lattice [13, HH . In many cases [3ll - l33l | , theoretical stud- 
ies have used the Bose-Hubbard model (BHM) or Time- 
Evolving Bloch Decimation (TEBD) approach, both of 
which become problematic when there are many atoms 
per well, as in the one-component experiments of fTTj | . 

For BECs in an optical lattice [34[ , phase diffusion and 
phase decoherence are fundamental quantum-mechanical 
phenomena [35U371 ] . In this paper, we study phase co- 
herence of interpenetrating two-component mixtures that 
are slowly loaded into relatively shallow state-dependent 
lattices, V _A(z,t) — sa^Er cos 2 (kz) for component (A) 
and V ,b = (or V 0j B(z,i) = s b (t) En sin 2 (kz)) for com- 
ponent (B) (si is a scale of lattice height for the compo- 
nent i and En is a recoil energy) , such that a single com- 
ponent would remain phase coherent after application 
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of the optical lattice. We construct a TWA simulation 
model for two-component BEC clouds which are indepen- 
dently phase-coherent in the initial state. Quantum and 
thermal fluctuations centered around mean-fields that in- 
duce the dephasing are represented by Bogoliubov quasi- 
particles in thermal equilibrium with a condensate gas. 
We focus on the effects of a second component which does 
not see an optical lattice directly (V ,b = 0) or alterna- 
tively undergoes a half-period mismatch optical lattice 
(Vo,B( z ,t) — SB(t)E R sm 2 (kz)). As the first case, we 
consider adding a second component to an original com- 
ponent of BEC, where an optical lattice that is effective 
for only the first component, V 0> a(z, t) is ramped up. The 
other component then sees an effective lattice from the 
mean-field potential associated with the first component. 
We find, as in the experimental studies mentioned above 
[l9l |20| , that the second component diminishes the phase 
coherence of the first component, and also experiences de- 
coherence itself relative to the initial fully coherent state. 
Because conditions differ (fD vs. 3D optical lattices, for 
example), we do not claim to have quantitatively valid 
models of these experiments, but instead we believe our 
results might serve as prototypical examples. 

In addition to the case of phase coherence loss over a 
range of density ratios of two components, we also con- 
sider a range of lattice heights (sa) for the component 
exposed to the optical lattice. We show that phase co- 
herence loss of both components is enhanced as the lattice 
height of one component is increased. 

Experimental studies on phase coherence properties of 
BECs have relied mainly on time-of-flight images taken 
after the release of traps and the analysis of their visi- 
bilities. The first-order correlations between two neigh- 
boring sites calculated here closely mimic the visibilities 
of the interference patterns of condensates released from 
optical lattices. 

For a qualitative explanation of the fragmentation 
processes described above, we adopt a simple Gaus- 
sian variational ansatz for single-particle Wannier func- 
tions. We find that the model shows a good agree- 
ment with the trend of fragmentation inferred from the 
above TWA calculations. To check further the va- 
lidity of this model, we investigate a different regime 
where the second component previously not exposed 
to a lattice is now fragmented by its own optical lat- 
tice, V ^B(z,t) — s b (t)Efi. sin 2 (kz) which is displaced 
a half-period from the lattice for the first component 
iy o ,A{z, t) — s a (t) Er cos 2 (kz)). Phase coherence loss 
of the first component under a fixed final lattice height 
will be shown to have non-monotonic dependence on frag- 
mentation of the second component, which is subject to 
variable lattice heights. 

The layout of this paper is as follows. In Sec. HB we 
construct the TWA model for ID two-component BEC 
beginning from a second-quantized effective Hamiltonian. 
In Sec. IIII| the TWA representation is applied to ini- 
tial states, where the Wigner probability distribution for 
the initial state is found, and we prepare an ensemble 



of initial states under the Wigner distribution. Sec. IIVI 
discusses relevant experimental parameter regimes and 
numerical methods for our calculations. We present the 
main results of the paper in Sec. [V] Sec. IV Al introduces 
a single state-selective optical lattice and shows the effect 
of an added component on phase coherence loss over a 
range of populational fractions of each component. We 
also implement a variational ansatz calculation to explain 
the patterns found above. Sec. IV Bl then continues the 
similar setup but with a fixed density and variable lattice 
heights to see the fragmentation induced by lattice height 
increase. In Sec. IV CI we find non-monotonic dependence 
of fragmentation on lattice height as two state-dependent 
peak-mismatched optical lattices are turned on. Finally, 
Sec. IVII is devoted to a few concluding remarks. 

II. DYNAMICS OF TWO-COMPONENT BECS 
IN THE TWA 

We consider a mixture of two Bose-Einstein conden- 
sates which is confined in a harmonic trap, where the 
two components are two hyperfine states of the same 
species [38] . The harmonic trap potential is Vh ar (x) = 
\m(uj 2 z 2 + uj 2 (x 2 + y 2 )) with a weak longitudinal trap 
frequency (u> z ) and a stronger transverse trap frequency 
(uj p ) (u> z < ui p ), so that the BECs are cigar-shaped. 

Assuming effective ID BECs with negligible transverse 
excitations, as explained in Sec. IIVI the effective ID two- 
component second-quantized Hamiltonian for the system 
is 

H= I dz$\(z)LSi(z) 

i=A,B •* 

+ 9u I dz$l(z)i>l(z)i>i(z)i>i(z) 

L i=A,B J 
+9AB J dz$ A (z)$ B (z)i)A(z)i>B(z), (1) 

and the Li is defined as, 
h 2 V 2 

L t = -- +V h ,i(z) + V ,i(z,t)- in. (2) 

Zrrii 

Here, we label the first species as 'A' and the second 
one as 'B'. For each species, m,i is the particle mass, /Lt« is 
the chemical potential, Vh,i(z) — miU! 2 z 2 /2 is the exter- 
nal harmonic trap potential, V ^(z, t) is the time- varying 
state-dependent optical lattice potential along the axial 
direction. For a 3D BEC, 

4 3D) = 2nh 2 aij (— + —), %,j = A, B; (3) 

where the g\f are intraspecies interaction strengths, 

9ij D ^' * 7^ J i s the interspecies interaction strength be- 
tween species i and j. an is the scattering length for the 
component H\ and a^; i ^ j is the interspecies scatter- 
ing length. For an effective ID BEC with a Gaussian 
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profile along the transverse direction, g^j — 2ftw p aij, if 
the masses are all equal. 

The equation of motion for the component i field, 

4>i(z), is 

. d 



dt 



•ipi = Hiipi = Litpi + ^2 9ij'4 , j'>Pji>i- (4) 



We construct a TWA method for the dynamics of two- 
component BECs under the nonequilibrium ramp-up of 
state-dependent optical lattices. We begin by project- 
ing the above two-component state onto its phase space 
within the Wigner representation. We thereby obtain 
a quasi-probability distribution function over the phase 
space, which is formulated to be analogous to the den- 
sity matrix in quantum mechanics [39| . Even though a 
positive-P representation is sometimes used for simula- 
tions of a quantum system, it is subject to instabilities, 
for example in highly populated modes Instead, we 
approximate the dynamics of one-dimensional trapped 
multimode BECs by considering the Fokker-Planck equa- 
tion in the truncated Wigner representation |40l |. 

The Fokker-Planck equation for the Wigner quasiprob- 
ability distribution yields the time evolution equation: 



dt 



i,j=A,B 



{(Li 



-9ij{ IV^O)l 2 - <kj)}i>i(z) 



9*3 



4 5ipi(z) S^j(z) 8ip*{z) 



(V>a,V>b) T , and d i3 - 



W(ip,ip*) + h.c, (5) 



1 (or 1/2) if i 



where ip = 
j (or i^j). 

The exact Fokker-Plank equation with the presence of 
the third-order term within the Wigner representation 
is difficult to solve both analytically and numerically in 
stochastic simulations (l6j . Therefore, the truncation in 
the TWA neglects the third-order derivative terms in Eq. 
([S]), which are smaller than the Gross-Pitaevskii first- 
order term in the total number, N. The second-order dif- 
fusion process term of usual stochastic processes, which 
can have a prominent role in enhancing fluctuations, is 
absent in the TWA. Also, in the TWA, the expectation 
values of those classical fields, Ylijtp*(zi)ipj(zj) in the 
Wigner representation correspond to the expectation val- 
ues of quantum operators that are symmetrically ordered, 
where the Wigner quasiprobability distribution function, 
W(ip, ip*), is a classical projection function corresponding 
to the density operator for the field operators in quantum 
mechanics: 

(Y[$i( z i)i>j( z j))w = 

i,3 

[ #$W$,$*)]lil>t(zi)M*i)- ( 6 ) 

J i,3 

We can translate the Fokker-Planck equation into the 
stochastic differential equation for the classical Wigner 



fields, ipi(z,t) |4ij . The resulting equation for a single 
realization of the Wigner fields that describes a single 
trajectory in phase space is 



ih 



dipj(z,t) 
dt 



[Li + ^2gij{\ipj{z,t)\ 2 - d i: 

3 

ipi{z,t), 



(7) 



Since the third-order diffusion process is neglected, 
the stochastic fluctuations during the time evolution are 
absent, but the initial state still has quantum fluctua- 
tions following the probability distribution given by the 
Wigner representation. Therefore, given the initial con- 
dition for each realization following the Wigner function, 
the classical field, ipi(z, t), evolves under the above deter- 
ministic trajectory which resembles the Gross-Pitaevskii 
equation except the small depletion terms indicated by 
the ll —dij" quantities in Eq. ([7}. 

The model we consider is assumed to have no loss 
mechanism from either two-body or three-body collisions. 



III. STOCHASTIC INITIAL STATES AND 
PHASE COHERENCE BETWEEN SITES 

For BECs with a large number of atoms at sufficiently 
low temperatures (T -C T c ), the Bogoliubov quasiparti- 
cle description (42j is a good approximation to the exact 
many-body dynamics of the system, provided that the 
number of noncondensate particles (N ex ) is sufficiently 
smaller than that of condensate atoms (N c ) (N ex <C N c ) 
43]. A more exact number-conserving theory would be 
based on an expansion in powers of i/Vn pi 0] 
usin g th e Particle-Number Conserving formalism (PNC) 
(46l |47| |. However, for N >> 1, the extra correla- 
tion terms between condensate and noncondensate modes 
would have a smaller effect than the Bogoliubov quasi- 
p article modes. 

In the Bogoliubov theory, the matter-wave field opera- 
tor, in addition to the condensate field operator, includes 
small quasiparticle amplitudes, 

tpi(z) = i/>io(z)&io + '^2(u ifl (z)a^ - VinW&l), (8) 

where d^o is the annihilation operator for the compo- 
nent i condensate mode, whereas d M is the quasipar- 
ticle annihilation operator for the collective mode fi. 
These operators satisfy the bosonic commutation rela- 
tion, [d^jdj,] = (5 p „, etc. The normalization conditions 
for the single-particle condensate amplitudes and for the 
Bogoliubov quasiparticle mode amplitudes are 



dzip* AO (z)ipAo(z) = J dziJj* BO (z)ipBo(z) = 1, (9) 

J dz[u* All (z)u A u{z)+U* BlJi {z)UBu{z) 
- v *An(z)VAv{z) - V* Bfl (z)VBu{z)) = 5^. (10) 
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The expectation values of the number operator corre- component and in the collective modes, 
spond to the populations in the condensate mode for each 

(or i0 a i0 ) = N i0 , (11) 

( ^ >= ^ = cxp( £ ,/It)-1< (12) 

where the Bogoliubov quasiparticles are in thermal equi- 
librium at T. 

The quasiparticle mode amplitudes satisfy the coupled 
two-component Bogoliubov-de Gennes equation: 



H A + h AA \ipA\ 2 

hAB^A^B 

-h A A(r A ) 2 

-hAB^A^B 



hABlpAlpB 
H B + h B B\i>B\ 

-hABlpA^B 

-h B B(r B ) 2 



-h-AAlpA 
-hAB^AlpB 

Ha + Haa^a? 

hABlpA'IpB 



-hABlpA^B 
-JlBBlpB 
tlABlpA^B 



>bb\^b 



UA^ \ 




1 e^UAp, 


UBf> 


- 


tfiUB^ 






-^VAp. 


VB^ ) 




\ -e^Bn 



(13) 



where /;, , = «/, , ^/N~N~ r 

As noted above, we generate a set of classical stochastic 
fields for the initial state sampled from the correspond- 
ing Wigner distribution function and obtain the dynam- 
ics of the system by averaging the statistical ensemble 
over individual trajectories in phase space. The expecta- 
tion values of symmetrically-ordered operators are calcu- 
lated by the weighted average of the corresponding clas- 
sical fields (ipw) with the Wigner distribution function, 
W(i^w, i^w) without extra modification terms. 

The classical stochastic fields are obtained by 
randomly-generated c-numbers via a given distribu- 
tion, a Ao, aso, c^m corresponding to quantum operators, 
&Ao,ctBo, and d M . Since the condensate mode opera- 
tors and quasiparticle mode operators commute and the 
component A and B condensate operators also commute, 
we independently sample the c-numbers. For the con- 
densate mode, the initial two-component superfluid state 
is approximated as a coexisting mixture of independent 
Glauber coherent states, where each coherent state pre- 
serves its own phase coherence. The Wigner function for 
the initial state is 



W(a,a*) = W A o(c(ao, u*Ao) w Bo(aBo, a B o) 
W B G(a,a*), 



(14) 



where a = (a^o, "bo, ol^) t and Wao , W B o , W B g are the 
Wigner distribution for the component A, the component 
B, and the Bogoliubov modes. The condensate mode 
Wigner distributions are given by (i = A, B) 



W i0 (ai ,a* ) = -exp 

7T 



2|a i0 



(15) 



where the ensemble averages are (ctio}w = VNio and 
(a* ctio)w = + i and W denotes the ensemble aver- 
age in the Wigner distribution. The distribution func- 
tion of coherent states has a Gaussian profile in the 
complex phase space with a variance of 1/2. For a 



large number of atoms (N 3> 1), quantum fluctuations 
around the mean classical field are relatively small, since 
AN/(N) = l/y/(N). Thus we can think of the initial 
state as a classical field with a small fluctuation in phase 
space. 

While the condensate modes have nonzero expectation 
values for the atom number of components, Ojo, the non- 
condensate mode operators, have zero expected popu- 
lations, for which the Wigner distribution is the product 
of uncorrelated Wigner functions for each mode [43j : 



W BG (a, a*) = JJW M (a^ap, 



W^(a M ,a*) 



ii 
2 



tanh 



exp 



2\a u \ tanh 



\knTJ 



(16) 



where Wbg(&, a*), W fl (a fl , a*) is the Wigner function 
for the total Bogoliubov modes and for each mode, re- 
spectively. The ensemble averages of quasiparticle modes 
satisfy the condition that they have the zero mean values 
and Gaussian variances, which broaden as the tempera- 
ture increases. 



(a*a v )w = S^[n u + -]. 



(17) 
(18) 



After the initial state of ID condensates is prepared 
in a harmonic trap, a state-selective optical lattice is ap- 
plied to only the A component along the axial direction. 
This will eventually lead to phase decoherence. We are 
especially interested in the short-range non-local coher- 
ence of subcondensates between neighboring sites. We 
define a subcondensate projection operator for each site 
I as in G3: 



/ th sitc 



dzipGp(z,t)tpi(z,t), 



(19) 
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where an is the annihilation operator for component i in 
the I th well, V'gp the solution of the GPE, normalized 
to one. The definitions of site positions are different for 
the two components as explained below. This operator 
is defined as a stochastic field operator whose amplitudes 
are projected over the ground state of each condensate 
mode. The projection method allows us to avoid compli- 
cated calculations of symmetrically-ordered multimode 

fields g. 

In this study, a state-dependent optical lattice for the 
component i is a sinusoidal function with a local max- 
imum at the center, V ,A{z,t) — s a (t) E r cos 2 (kz) (or 
V ,B{ z ,t) = s b (t) E R sin 2 (kz)), where is the scale of 
lattice height for the component i, and En = h 2 k 2 /2m 
is the recoil energy. For the lattice with a wavelength A 
(k = 2tt/A), the periods of the lattices are all equal to 
d = A/2. There is no optical lattice for component B for 
most of the calculations in this study. 

The ground state at each instance of time t is numer- 
ically obtained by a standard imaginary time evolution. 
The stationary ground state wavefunctions of both com- 
ponents are symmetric with respect to the origin, i^a 
and tpB are localized at the odd sites z — ±(2n + l)d/2, 
at the even sites z = ±2nd/2 (n = 0,1,2,...), respec- 
tively. Repulsive interspecies interactions repel compo- 
nent B atoms from the localization sites of the component 
A. Also, we selectively define the projection operators of 
the component A and B only for the odd and even sites 
in Eq. rjj 

We now consider moments of the Wigner function of 
interest. First, the occupation number of component i in 
the I th site is 

Tin = (afjflti) = {dl^ijw ~ \ (20) 

For a BEC with a large number of atoms, the trans- 
formation offset 1/2 between the symmetrically-ordered 
operators and the normally-ordered operators can be dis- 
regarded. 

A more interesting observable is the degree of coher- 
ence. The general n th order correlation function is de- 
fined as 

(„) = (a\{h)al(t2)---aUtn)an(t n )...a2{t2)ai(ti)) ^ 

nr<«!(^(a i/2 

where a, denotes any projection operator. The equal- 
time first-order coherence is the phase coherence of com- 
ponent i between two sites at the time t: 

(i) = c (t)= l(4(*)a«'(*))l 

* (<4(*)M*))<4'(*)M*))) 1/2 

= iM, (22 ) 

where in the last equation the notation is simplified via 
Eq. [501 For brevity, we now omit the time dependence 
from the expectation values of the condensate mode op- 
erators. 



IV. NUMERICAL METHOD 

The system we discuss is a two-component ID BEC 
confined by the same harmonic trap. In this discus- 
sion, we consider the two hyperfine states of 87 Rb atoms, 
\F = l,m F = -1) and |F = 2,m F = -2). Since the dif- 
ferences between interaction strengths (aAA,o-BBi a AB) 
are small, we assume that the two atoms share the 
same intraspecies and interspecies interaction strength 
(a s = ciaa — &bb — a AB = 5.5nm) and also assume 
that they have the same masses (tuaa = ttibb = m). 

The condensation of cigar-shaped ID atomic clouds is 
achieved in an anisotropic harmonic trap with the tight 
confinement along the transverse direction (uj z = 2tt x 
l30Hz,uj p = 2ir x 2.71kHz) and the aspect ratio is 21. 
In this work, the number of atoms ranges from 0.5 x 10 3 
to 5 x 10 3 . The optical lattices are generated by red- 
detuned off-resonant lasers with a wavelength A = 785 
nm, so that the period is d = A/2 = 392.5 nm, and the 
recoil frequency is u>r = (ft/2m)(27r/A) 2 = 27rx3.73 kHz. 

The one-dimensionality of cigar-shaped BECs in the 
harmonic trap at T = is achieved when l z > 
£ > l p or fan < hjjp [48[, where l z — (h/mujz) 1 / 2 
and l p = (h/mojp) 1 / 2 are the longitudinal and the 
transverse zero-point oscillation length respectively, and 
£ = (l/47rri3Da s ) 1 ' 2 is the healing length, and Hsn = 
h 2 /2m(l5Na s /l 2 lp) is the chemical potential correspond- 
ing to the interaction energy. Furthermore at T > 0, l p 
is required to be smaller than the de Broglie wavelength 
At {lp < A T ), where A T = {2nh 2 /mT) 1 / 2 For the 

87 Rb atoms in a trap with frequencies given above, the 
ratios, Z p /£ < 2 and hzd/Tu^p < 3. 

Each BEC component lies in the Thomas- Fermi regime 
between the full 3D dynamics and the true ID dynam- 
ics with transverse excitations almost frozen out. Even 
though the BEC is in the crossover between 3D and ID, 
the low-energy excitation modes in 3D are effectively ID 
provided that the temperature is sufficiently below the 
energy of the transverse oscillator (T < hco p ) which is 
the case here (T = 0) 50]. The Thomas- Fermi radius 
ranges from 5 ~ 121 z . 

The dimensionless coupling strength of interaction en- 
ergies in this work is 7 = mgm/h 2 niD < 2 x 10~ 3 
and the reduced temperature is r = 2mkBT /h 2 n\ D = 
[5lj ]. Therefore, the ID Bose gas can be effectively de- 
scribed by the Gross-Pitaevskii equation in the regime 
( r2 ^ 7 ^ 1) far from the Tonks-Girardeau regime 
(7 ^ !)• The nonlinearity gmN /fmj z l z [HI ranges from 
120 to 1200. 

The numerical preparation of initial states requires 
ground state wavefunctions, the Bogoliubov quasipar- 
ticle excited modes, and their stochastic distributions 
governed by the Wigner functions. We find the ground 
state wavefunctions by numerically integrating the GPE 
in imaginary time with a time step of coftSt = 0.005 with 
3072 spatial grid points along the axial direction. We 
utilize the second-order split-operator method to inte- 
grate the time evolution of wavefunctions, in nonlinear 
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as well as linear regimes. Using the ground state solu- 
tions of the GPE, we obtain quasiparticle wavefunctions, 
UAn{z),UBn{z),VA t i{z),VB li {z) for energy, by the di- 
agonalization of the Bogoliubov-de Gennes equation [Eq. 

ES. 

We calculate the ensemble average of stochastic fields 
along the trajectory and find their coherence. Stochastic 
quantum fluctuations are appended to the initial mean- 
field state for generation of the ensemble of Wigner- 
distributed initial states, in which step we perform the 
Gaussian random variable generation of order parame- 
ters (aio,aj) for both components. For the condensate 
mod e, th e mean of ai is v^io and its width of deviation 
is y/T/2, whereas for the Bogoliubov quasiparticle mode, 
the mean of is zero and the width is ^/l/2 for T = 0. 
A single sample of stochastic fields, 4>w{ x ) is obtained by 
configuring the wavefunction profiles with the generated 
stochastic order parameters. 

The condition for numerical validity of the TWA 
method in the Bogoliubov theory is that the condensate 
mode must be highly populated compared to the noncon- 
densate mode so that the quantum fluctuation is small 
being dominated by the condensate field. In other words, 
the TWA in the mean-field theory is valid with a rela- 
tively small number of excited modes compared to the 
number of condensate particles in the system, N ^> M/2, 
where N is the total number of atoms, M is the number of 
Bogoliubov quasiparticles. This is a regime different from 
other exact numerical methods, for example, the Time 
Evolving Block Decimation (TEBD) method or the Den- 
sity Matrix Renormalization Group (DMRG) with the 
Bose-Hubbard Model, in which cases each site is limited 
to a low filling factor since the Hilbert space increases 
exponentially with the number of atoms and the number 
of sites. 

We perform the simulation with an ensemble of states 
consisting of 500 samples for the TWA distribution func- 
tion to achieve sufficient convergence. The time evo- 
lution of ensembles has the typical time step given by 
Lu R 5t = 0.005, i.e. St = 0.2 (is. 



V. PHASE DECOHERENCE AND 
FRAGMENTATION OF TWO-COMPONENT 
BECS 

A. Increasing fragmentation of A atoms, decreasing 
fragmentation of B atoms as the B fraction is raised. 

We now examine the phase decoherence patterns of 
two-component BECs at T = 0, driven by a single state- 
dependent optical lattice. Both components are trapped 
by the same anisotropic harmonic potential. Without 
loss of generality, component A in this work is placed in 
an optical lattice whereas no external lattice is applied 
to component B. In this subsection, we fix the total atom 
number and the ramp-up time, and vary the number ra- 
tio of A to B atoms. Both species fragment as the A 




co R t/10 2 co R t/10 2 

FIG. 1. Phase coherence of component A (on left) and 
component B (on right) for various populational fractions 
(0 < f B < 1), for N = 5 x 10 3 , t ru = 250/uj r = 11ms 
for 87 Rb. The numerical simulation shows the results from 
when lattice loading begins. The fractions of component B, 
/b, and the line colors are (in the large t limit, top to bottom 
on the left; bottom to top on the right): 0, black (for left 
side only); 0.1, red; 0.2, green; 0.4, blue; 0.6, pink, (online in 
color) 



lattice is applied, but with more B atoms, the B atom 
fragmentation is much reduced. 

Having prepared the initial state of superfluid BECs 
placed in the harmonic trap, we linearly turn on the op- 
tical lattice up to a final height of s maXt A = 10 m ramp-up 
time of uirTru — 250 (tr;j=11 ms), then maintain the 
height until the end of simulations: 

V a ,A(z,t) = s A {t)E R cos 2 (kz), 
V o , B (z,t)=0, (23) 

where s A (t) = s max , A t/T RU , < t < t ru . 

We fix the total number of atoms, N to t = 5 x 10 3 , 
and vary the fractions of components A and B, j a and 
Jb = 1 — f A , in order to see the effects of interspecies 
interaction and imbalanced populations on phase deco- 
herence. With this value of Ntot, an d the trap and lattice 
properties given above, we obtain as many as 75 atoms in 
the central well of the quasi-lD lattice. Experimentally, 
working also with 87 Rb, Campbell et al. [52j were able 
to put at most 5 atoms per site in their 3D optical lat- 
tice. Because of the tighter confinement in [52], 75 atoms 
per well in our simulations would actually correspond to 
about 8 atoms per well in [52[ for the same density at the 
peak. In actual experiments, the total number of atoms 
would need to be reduced over the value used here. This 
would lead to a reduction of the following coherence loss 
effects. 

From previous experimental and theoretical work (Til 
ITijj ]. we expect component A in the absence of B atoms 
to exhibit phase decoherence under certain conditions. 
As the periodic lattice rises into the BEC cloud, the re- 
gions occupied by the lattice peaks are locally avoided 
by ground state component A and the wavefunctions are 
eventually fragmented to some degree. The interaction 
energies increase and the tunneling rate of the wavefunc- 
tions between the adjacent sites is reduced so that the 
fluctuation in each subcondensate breaks the long-range 
phase coherence. 
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FIG. 2. The profile of combined potential for the component 
A, V c ,a(z) = Vh,a{z) + V ,a(z) with s max ,A = 10 and for 
the component B, V c ,b{z) = V hjB (z) + g ab\4> a{z)\ 2 . The 
population ratio (q) is 4 and Na = 4 x 10 3 , JVb = 1 x 10 3 . 



Figure [JJ shows the change in phase coherence between 
the center and nearest neighbor well, Ci-o±, for compo- 
nent A (left) and B (right) from the time the lattice be- 
gins to ramp up, to a large time limit. The coherence 
changes for other distant wells (C, ; o2 5 Ci ; 03, Ci ; o4, etc.) 
exhibit a similar pattern, but with more coherence loss 
at a given time. The first-order correlation functions be- 
tween sites are closely related to the visibilities of the 
interference pattern when the atoms are released under 
gravity [H, [53J. For one-component BEC cases, a com- 
plete loss of phase coherence would imply a transition to 
the Mott insulator state. In these calculations, the max- 
imum lattice height does not reach the Mott insulator 
regime, as indicated by the observation that in Fig. [TJ 
CU ; oi remains very close to unity if /a = 1. However, as 
/b increases, component A exhibits decoherence. 

A new feature in this two-component case is the reduc- 
tion of phase coherence of component B, which is induced 
as for component A but with the role of the optical lat- 
tice replaced by the atomic mean-field potential formed 
by component A's periodic localization. In the GPE for 
the component B, 

_dtp B (z,t) 



ih- 



dt 



[Lb + 9BB\^B{z,t)\ 2 

+gAB\*l>A(z,t)\ 2 ]lp B (z,t), 



(24) 



such spatial variation of potential is expressed by the 
term gAB\i>A{x, t)\ 2 . In Fig. [2j we show the optical 
lattice with the harmonic trap, which directly affects 
coherence properties of component A, and the mean- 
field lattice of A with the same harmonic trap, acting 
on component B, for the case N = 5 x 10 3 , Na — 
4 x 10 3 . The distortion by the harmonic trap poten- 
tial is almost negligible around the center. The depth 
of the optical lattice and the interaction strength of 
the mean-field lattice are comparable ((I a i,A{z)\max — 

Ial.A(z)\ m in) / (V ,A(z,t)\ max - V a ,A [z, *) \min) — 0-6) for 

}b — > as can be verified by Fig. [3J for which 



IaLA{z) = gAB\i>A{z)f 



(25) 



Due to the presence of the mean-field lattice, the tunnel- 
ing amplitude between the localization sites for compo- 
nent B is reduced, resulting in coherence loss, as shown 
on the right of Fig. [TJ 



The phase decoherence of component A is greater in 
the presence of component B than without component B, 
and increases as fs increases. Note also that the mean- 
field potential from B atoms acting on A atoms is in phase 
with the optical lattice, and thus effectively raises the pe- 
riodic potential that A atoms see, therefore contributing 
to the loss of coherence of the A atoms. However, com- 
parison with the degree of coherence for A atoms alone 
as a function of lattice height shown in the next section, 
elevation of the effective lattice acting on A atoms does 
not explain fully the decrease of coherence shown in Fig. 
[JJ (left side). Evidently the stochastic nature of the atom 
distributions also plays a role. 

The experiment in 20] has shown a similar dependency 
on populational fractions but with two peak-matched 
state-dependent optical lattices in order to place two 
components at the same lattice site. 

To gain another perspective on this process, we expand 
the wavefunctions in an array of Wannier-like orbitals, 

Wi(z), 



auWi(z - Ru), 



(26) 



where the single particle wavefunction, Wi(z-Ru) is cen- 
tered at R A i = (21 + l)d/2, R B i = 2ld/2 for each com- 
ponent, and the operators ay satisfy the bosonic com- 
mutation relation, [a^,aj ; ,] = Sij5u>. The variationally 
minimum solution of orbital wavefunctions implicitly de- 
pends on the occupation per site. Putting this set of 
orbitals into the Hamiltonian in Eq. (fT]), we obtain 



H = - Ji-w (a' u au' + d\ v aa) + Uadl^aaau 
+ UABa Al a Bl aAiaBi + eua^au, (27) 



id 



where 



JidV = — I dzwi(z - Rii) 

h 2 



2m,i 



V +V h ,i(z) + V ,i(z,t) 



U tj = g I dzw 2 (z - Ru)w 2 (z - Rji). 



Wi(z - Rii>), 
(28) 



In the tight-binding limit, the Wannier functions can 
be written as Gaussian functions p54j . When the tight- 
binding limits (V 0iA (z,t)\ max -V 0i A(z,t)\ min ) > E R and 

{hlAz)\ max ^al.A ruin 

) 3> Er are achieved, the high 
vibrational modes for each component are not occupied 
especially at the initial temperature T — so that the 
profile of each component can be well described by the 
ground state, the Gaussian wavefunction. Starting from 
the initial trial state of infinite ID BECs {V^h — 0) in the 
periodic state-dependent optical lattice, we employ the 
Gaussian variational ansatz for a single-particle orbital 
placed on each site, Wi(z — Ra) = (1/na 2 ) 1 ^ 4 exp(— (z — 
Ru) 2 /2a 2 ), with the density of atoms per site equal to 
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FIG. 3. The on-site interaction energies, Ua (solid line) and 
the widths of single-particle wavef unctions, Oi (dashed line) 
for component A (left) and for component B (right), show- 
ing opposite changes as the impurity component B populates 
increasingly. Here SA,max=W, N to t = 5x 10 3 . 
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FIG. 4. The effect of final optical lattice heights on phase 
coherence of component A (on left) and component B (on 
right). The final lattice height (sa) for each curve is following: 
the black line (3), the red line (6.5), the green line (10), the 
blue line (13.5), and the pink line (17). 



the average density of the center site calculated from the 
GPE (uu = n\ ) ■ Here, the widths of Gaussian wave- 
functions are variational parameters, as in [55H58l |. 

Within the Gaussian approximation, we obtain the 
minimized Gross-Pitaevskii energy functional, where the 
interaction energies and the tunneling amplitudes can be 
calculated from variational parameters. The local single- 
band Gaussian state is known to be accurate for the cal- 
culation of on-site interaction energies even for shallow 
lattices (~ 3-Er) with the overlap between the true Wan- 
nier function and the Gaussian wavefunction nearly equal 
to 1.0 [|{| • Since the Gaussian ansatz can be quite impre- 
cise for the calculation of tunneling amplitudes because 
of the tail of Gaussian functions |59l |60|, however, we 
concentrate on calculating on-site interaction energies. 

In Fig. |31 we show on-site interaction energies for com- 
ponent A and B as a function of the fraction of compo- 
nent B (/b) using the same parameters as in the TWA 
simulations, s max ,A = 10 and N to t = 5 x 10 3 . Higher 
on-site energies, Uu correspond to greater localization, 
(smaller a) and reduced nonlocal coherence, C^oi- 

The bosonic mixture with a high population ratio, 
q 3> 1 (/a — 1), can be approximated by the foreground 
component A with a B impurity. Within the first-order 
approximation, the average strength of interspecies inter- 
action (~ gABfAN) over its spatial variation, is greater 
than when j a < 1< The interaction strength varies over 
space because of the component A's modulational vari- 
ance and the points of localization sites for two com- 
ponents differ by d/2 because of repulsive interaction. 
Therefore, we can qualitatively interpret the strength- 
ened phase decoherence of component B for high q (large 
Na/Nb) by the increased strength of the mean-field lat- 
tice formed by the component A. 

On the other hand, the larger phase decoherence of 
component A as Ja — > can be understood by the com- 
bined effects of a larger number of B atoms and also of 
the broadening of component B distribution in each B 
well and the narrowing of the A distribution in each A 
well. Due to the repulsive nature of interspecies inter- 
action, the minimized Hamiltonian is found in the bal- 
ance between reducing the spatial overlap of the two- 
components' amplitudes and weakening the intraspecies 
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16 



0.75 
0.50^ 



^0.25 
20 



FIG. 5. The on-site interaction energies, Ua (solid line) and 
the widths of single-particle wavefunctions, <r; (dashed line) 
for component A (Uaa on left) and for component B (Ubb on 
right) for Ntot = 5 x 10 3 . Both components become more frag- 
mented as the state-dependent optical lattice for component 
A only has been amplified more. 



interaction energies of each component. The competition 
between the interaction strengths of two components gen- 
erates one component's localization and the other com- 
ponent's delocalization as the population ratio q changes. 



B. Increased fragmentation of both A and B atoms 
when the optical lattice height is raised 

In this section, we fix the fractions of the two compo- 
nents and show how the phase coherence changes as a 
function of time, depending on the final lattice height for 
component A. As in Sec. IV A[ V ,B(z,t) — 0, but now 
the atom numbers are fixed at Na = 4.0 x 10 3 ,A^b = 
1.0 x 10 3 . As before, the optical lattice for component 
A linearly increases up to the indicated value of s max ,A'- 
the ramping-up time is ujrtru = 100 in this case. 

The changes in phase coherence, Ci ; oi, are shown in 
Fig. @]for component A and B. As the lattice becomes 
deeper, the phase coherence of component A decreases 
as expected, since it is fragmented by the external lattice 
height increase even without consideration of interspecies 
effect. Decoherence of component B is enhanced as well 
because of the growth of the mean-field lattice from com- 
ponent A. In Fig. [SI the on-site interaction energies are 
displayed for both components in the regime where the 



two components are sufficiently localized (s r 



> 



10). 
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As shown in Fig. [3l an increase in on-site energy implies 
a smaller a, hence tighter localization within the effec- 
tive well. Thus Fig. [5] indicates, as expected, that the 
degree of localization is higher for deeper lattice heights. 
Comparing results between the two components in Fig. 
13 the localization of component A is evidently stronger 
than component B for s max .A > 10, which explains the 
greater phase decoherence in component A than in com- 
ponent B in Fig. |4j 

As is evident from comparing Figs. [5]and[3l when the 
A lattice height rises, the exchange of spatial occupation 
between the two components that has been observed in 
Sec. IV Al does not occur. Instead, component B's local- 
ization is strengthened as well as component A's. The 
atomic mean-field lattice formed by component A be- 
comes more peaked as it enters into the Mott insulator 
regime so that the tunneling amplitude of component B 
between its neighbouring sites is reduced. The loss in 
first-order spatial correlation between wells can be in- 
duced by increasing the height of barriers [6l| , which are 
provided by atomic mean-field lattices in this case. 



C. Non-monotonic fragmentation (A) - 
Fragmentation (B) 

Up to this point, component B has not been subjected 
directly to an optical lattice, but is localized simply by 
interaction with the mean held resulting from component 
A, and the interspecies interaction. Additional insight 
into the localization process can come from applying to 
component B an optical lattice that is out of phase with 
the lattice applied to component A so as to strengthen 
the localization effect on B atoms in addition to the mean 
held produced by component A, acting on component B. 

When an optical lattice is applied to only one com- 
ponent as in the previous sections, the decoherence as a 
function of lattice height has been monotonic. In this 
section, we examine the dependence of phase decoher- 
ence of one component on the fragmentation of the other 
component in the presence of two peak-mismatch optical 
lattices as in [62, 63]. 

To the BEC mixture with the asymmetric population 
ratio (q = 0.25, N A = 1.0 x 10 3 ,N B = 4.0 x 10 3 ), we 
gradually apply two state-dependent optical lattices: 



V 0tA (z,t) = s A (t)E R cos 2 (kz), 
V , B (z,t) = s B (t)E R sm 2 (kz), 



(29) 



where Sj(i) = s maXi it/Tmr- The final lattice height for 
the component A is s max ,A = 10 and the final lat- 
tice height for the component B is a variable parame- 
ter in different simulations ranging from s max .B = 3 to 
s m ax,B = 17 and all other conditions are the same as in 
the previous section. 

In Fig. |51 we show the change of phase coherence while 
varying the lattice height for component B. The drop in 
phase coherence of component B reflects the expectation 
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FIG. 6. The effect of localization induced by B on phase 
coherence on component A (on left) and component B (on 
right). The final lattice height for B (sb) is following: the red 
line (3), the green line (6.5), the blue line (10), the pink line 
(13.5), and the yellow line (17). 
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FIG. 7. The on-site interaction energies, Uu (solid line) and 
the widths of single-particle wavefunctions, a (dashed line) 
for component A (Uaa on left) and for component B (Ubb 
on right) as a function of the lattice height for B (sb). 



that higher optical lattices induce more coherence loss be- 
tween two neighboring sites, neglecting interspecies inter- 
actions for the moment. And also as expected, a higher 
lattice for component B leads to larger on-site interac- 
tion energies, as shown in Fig. and this in turn implies 
more localization for component B, for reasons discussed 
above. 

The new phenomenon here is the dependence of phase 
coherence of component A on optical lattice heights for 
the component B. As can be seen from Fig. Ufa), the 
phase coherence of component A at the maximum time 
shown is diminished as s max .B increases from 3 to 10, 
but then rises again for s ma x,B increases beyond 10. This 
non-monotonic dependence on the other component's lat- 
tice height can also be seen in the on-site energy plotted 
in the left panel of Fig. [7j which shows a maximum of 
Uaa at s ma x,B ~ 15, and then a small decrease. In ad- 
dition to the expected fragmentation due to the optical 
lattice applied directly to component A, the deepening 
mean-field potential from component B acting on com- 
ponent A further reduces the tunneling of component A 
between adjacent sites until the tunneling of component 
A suppressed by increased mean-field lattice height of 
component B becomes gradually freed by decreased wall 
width of the mean- field lattice (erg). 
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VI. CONCLUSION 

In summary, we have investigated nonlocal intraspecies 
phase coherence of two-component BECs under the grad- 
ual loading of state-dependent optical lattices. We have 
used the TWA to model the dynamical behavior of two 
species of atoms in periodic potentials, in particular to 
calculate the reduction of phase coherence between atoms 
in the central and nearest neighbor site, C^oi, i = A or B. 
The initial state is a stable mixture of overlapping two- 
component BECs in a harmonic trap. An optical lattice 
is applied to the A atoms but not directly to the B atoms. 
Both species are fragmented into localized wells differing 
by a half-period between the two components. The B 
atoms are fragmented by the effective barriers produced 
by interaction with the spatially varying distribution of 
A atoms. 

First, with a constant total atom number and constant 
lattice amplitude, we have observed that when the frac- 
tion of B atoms increases, then in the long-time limit, the 
fragmentation of A atoms increases, but the fragmenta- 
tion of B atoms decreases, consistent with experimental 
observations in [l9|, [2(j. The increasing fragmentation 
of A atoms is associated with higher effective barriers 
produced by accumulation of B atoms spatially out of 
phase with the A atoms. Alternatively, the Gaussian 
variational wavefunction ansatz provides additional in- 
sight into the fragmentation mechanisms. Since intra- 
and inter-species interactions are all repulsive, when the 
wavefunction of the B species broadens within each lo- 
calization site due to increasing numbers of B atoms, the 
B atom fragmentation is reduced, and the increased den- 
sity repels the A atoms from the high-density regions of 
the B atoms, thereby enhancing the localization of the A 
atoms and diminishing site-to-site phase coherence. 

Secondly, we have seen increasing fragmentation of 
both A and B atoms with higher lattice heights for A. 



In this case, an external optical lattice drives localiza- 
tion and phase decoherence of A atomic fields which are 
directly exposed to lattice wells, whose mean- field poten- 
tial in turn forms an effective lattice for B atoms, reduc- 
ing its phase coherence between sites. We conclude that 
localization effects in deeper lattices play a key role for 
both components in reducing phase coherence. 

Finally, when a lattice is applied also to the B atoms, 
then as a function of its amplitude, non-monotonic frag- 
mentation of A atoms occurs; the A atom on-site in- 
teraction energies reach a maximum, and then decrease 
when the B atoms become more and more localized. This 
shows in a dramatic way how the effect of increased lo- 
calization (or equivalently, on-site energy) of one species 
acting on the other, can actually saturate. 

All the calculations pertain to a situation in which 
the two species are different hyperfine states, with equal 
masses and equal inter- and intra-species scattering 
lengths. Since this work is exploratory in nature, we have 
not attempted an extensive survey of parameter space by 
varying atom numbers, masses, and scattering lengths. 
We suggest that effects similar to what we obtain could 
occur if the two species were actually different types of 
atoms (Rb and K, for example), such as has been ob- 
tained in recent experiments. 
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